Instalamos los paquetes necesarios
library(dplyr)
library(sf)
Cargamos los datasets
estaciones <- read.csv("data_smn/preprocessed/estaciones_smn_v2.csv")
horarios <- read.csv("data_smn/preprocessed/datohorario20210420.csv")
Observamos como esta compuesto el dataset
glimpse(estaciones)
Rows: 123
Columns: 9
$ NOMBRE <chr> "BASE BELGRANO II", "BASE CARLINI (EX JUBANY)", "BASE ESPERANZA", "BASE MARAMBIO", "BASE ORCADAS", "BASE SAN MARTI…
$ PROVINCIA <chr> "ANTARTIDA", "ANTARTIDA", "ANTARTIDA", "ANTARTIDA", "ANTARTIDA", "ANTARTIDA", "BUENOS AIRES", "BUENOS AIRES", "BUE…
$ LATITUD_GRADOS <int> 77, 62, 63, 64, 60, 68, 36, 38, 37, 36, 34, 38, 37, 36, 34, 34, 34, 34, 36, 37, 34, 34, 34, 35, 36, 35, 37, 35, 34…
$ LATITUD_MINUTOS <int> 52, 14, 24, 14, 45, 8, 50, 44, 43, 12, 32, 1, 26, 21, 36, 49, 33, 58, 2, 56, 33, 41, 40, 27, 53, 52, 36, 22, 27, 1…
$ LONGITUD_GRADOS <int> 34, 58, 56, 56, 44, 67, 59, 62, 59, 61, 58, 61, 61, 57, 58, 58, 60, 57, 59, 57, 58, 58, 58, 60, 60, 61, 62, 57, 58…
$ LONGITUD_MINUTOS <int> 34, 38, 59, 43, 43, 8, 53, 10, 47, 4, 40, 20, 53, 44, 36, 32, 55, 54, 8, 35, 49, 44, 38, 53, 13, 54, 23, 17, 35, 1…
$ ALTURA <int> 256, 11, 24, 198, 12, 7, 147, 83, 207, 94, 26, 247, 233, 9, 12, 20, 81, 23, 36, 21, 32, 36, 24, 76, 166, 87, 304, …
$ NRO <int> 89034, 89053, 88963, 89055, 88968, 89066, 87641, 87750, 87649, 87640, 87570, 87683, 87637, 87648, 87571, 87576, 87…
$ NroOACI <chr> "SAYB", "SAYJ", "SAYE", "SAWB", "SAYO", "SAYS", "SAZA", "SAZB", "SAZJ", "SAZI", "SADO", "SABP", "SAZC", "SAZD", "S…
glimpse(horarios)
Rows: 2,041
Columns: 8
$ FECHA <int> 20042021, 20042021, 20042021, 20042021, 20042021, 20042021, 20042021, 20042021, 20042021, 20042021, 20042021, 20042021, 2004…
$ HORA <int> 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, …
$ TEMP <dbl> 20.7, 19.8, 19.5, 19.2, 19.2, 18.9, 19.0, 18.8, 19.0, 19.1, 19.9, 20.3, 21.8, 22.9, 23.2, 23.7, 23.2, 23.4, 22.5, 22.5, 22.1…
$ HUM <int> 70, 79, 79, 81, 81, 84, 84, 84, 84, 84, 81, 82, 68, 72, 68, 66, 71, 68, 75, 72, 73, 68, 68, 68, 90, 88, 91, 90, 93, 94, 98, …
$ PNM <dbl> 1020.8, 1020.9, 1020.9, 1020.4, 1020.0, 1020.3, 1020.7, 1021.1, 1021.6, 1022.2, 1022.4, 1021.9, 1021.7, 1021.0, 1020.1, 1019…
$ DD <int> 70, 70, 70, 70, 70, 50, 70, 50, 50, 50, 50, 50, 50, 70, 70, 70, 70, 70, 90, 90, 90, 90, 70, 70, 50, 50, 50, 20, 50, 20, 50, …
$ FF <chr> "28", "28", "26", "26", "22", "22", "17", "17", "15", "19", "20", "24", "20", "19", "20", "17", "15", "15", "19", "20", "22"…
$ NOMBRE <chr> "AEROPARQUE AERO", "AEROPARQUE AERO", "AEROPARQUE AERO", "AEROPARQUE AERO", "AEROPARQUE AERO", "AEROPARQUE AERO", "AEROPARQU…
Observamos que la variable que representa la velocidad del viento en km/h es de tipo Char, por lo cual debemos convertila en int.
horarios$FF <- as.numeric(horarios$FF)
NAs introduced by coercion
summary(estaciones$ALTURA)
Min. 1st Qu. Median Mean 3rd Qu. Max.
3.0 53.5 145.0 331.5 455.0 3459.0
summary(horarios$HORA)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.00 7.00 12.00 12.16 18.00 23.00
summary(horarios$TEMP)
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
-20.00 16.10 20.20 19.09 23.80 33.60 1
summary(horarios$HUM)
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
18.00 59.00 72.00 71.18 85.00 100.00 4
summary(horarios$PNM)
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
960.8 1011.0 1015.1 1075.3 1019.2 3133.0 23
summary(horarios$DD)
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.0 50.0 70.0 118.9 160.0 990.0 1
summary(horarios$FF)
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.00 9.00 15.00 16.06 22.00 107.00 1
hist(estaciones$ALTURA, main = "Histograma de la altura", xlab = "Altura")
hist(horarios$HORA, main = "Histograma de horarios", xlab = "Horarios")
hist(horarios$TEMP, main = "Histograma de temperatura", xlab = "Temperatura")
hist(horarios$HUM, main = "Histograma de humedad", xlab = "Humedad")
hist(horarios$PNM, main = "Histograma de presion atmosferica", xlab = "Presion Atmosferica")
hist(horarios$DD, main = "Histograma de la direccion del viento", xlab = "Direccion del viento")
hist(horarios$FF, main = "Histograma de la velocidad del viento", xlab = "Velocidad del viento")
Creamos unos boxplot para visualizar la distribucionde datos que potencialmente nos interecen para proceder con el analisis
boxplot(estaciones$ALTURA, main = "Boxplot de la altura", xlab = "Altura")
boxplot(horarios$HORA, main = "Boxplot de horarios", xlab = "Horarios")
boxplot(horarios$TEMP, main = "Boxplot de temperatura", xlab = "Temperatura")
boxplot(horarios$HUM, main = "Boxplot de humedad", xlab = "Humedad")
boxplot(horarios$PNM, main = "Boxplot de presion atmosferica", xlab = "Presion Atmosferica")
boxplot(horarios$DD, main = "Boxplot de la direccion del viento", xlab = "Direccion del viento")
boxplot(horarios$FF, main = "Boxplot de la velocidad del viento", xlab = "Velocidad del viento")
library(DataExplorer)
plot_boxplot(horarios, by = "HORA")
plot_boxplot(horarios, by = "TEMP")
Ahora, veamos que tan correlacionadas estan estas variables
library(ggcorrplot)
corr <- cor(horarios[, c(2,3,4,5,6,7)], use = "complete.obs")
ggcorrplot(corr, type = "lower", outline.col = "black",
lab=TRUE,
ggtheme = ggplot2::theme_gray,
colors = c("#6D9EC1", "white", "#E46726"))
Como correlacionan estas variables con la altura?
Hay datos nulos en estos datasets?
print("Nulos en horarios")
[1] "Nulos en horarios"
which(is.na(horarios))
[1] 4853 6710 6894 7999 8000 8401 8935 9137 9138 9139 9140 9141 9142 9143 9144 9556 9557 9558 9559 9560 9561 9562
[23] 9563 9564 9565 9566 9567 9605 10976 13017
print("Nulos en estaciones")
[1] "Nulos en estaciones"
which(is.na(estaciones))
integer(0)
Borramos los nulos del dataset
estaciones = na.omit(estaciones)
horarios = na.omit(horarios)
# Analisis de simetria
library(psych)
skew(horarios$FF)
[1] 1.591713
kurtosi(horarios$FF)
[1] 7.23501
Necesitamos convertir las latitudes y longitudes a
library(biogeo)
latitud <- c()
longitud <- c()
for(i in 1:nrow(estaciones)) {
latitud[i] <- dms2dd(dd = estaciones[i, "LATITUD_GRADOS"], mm = estaciones[i, "LATITUD_MINUTOS"], ss = 0, ns = "S")
longitud[i] <- dms2dd(dd = estaciones[i, "LONGITUD_GRADOS"], mm = estaciones[i, "LONGITUD_MINUTOS"], ss = 0, ns = "S")
}
estaciones['LATITUD'] <- latitud
estaciones['LONGITUD'] <- longitud
Unimos las dos tablas usando la variable NOMBRE como punto para combinar los datasets
data <- full_join(estaciones, horarios, by = c("NOMBRE" = "NOMBRE"))
glimpse(data)
Rows: 2,040
Columns: 18
$ NOMBRE <chr> "BASE BELGRANO II", "BASE BELGRANO II", "BASE BELGRANO II", "BASE BELGRANO II", "BASE BELGRANO II", "BASE BELGRANO…
$ PROVINCIA <chr> "ANTARTIDA", "ANTARTIDA", "ANTARTIDA", "ANTARTIDA", "ANTARTIDA", "ANTARTIDA", "ANTARTIDA", "ANTARTIDA", "ANTARTIDA…
$ LATITUD_GRADOS <int> 77, 77, 77, 77, 77, 77, 77, 77, 62, 62, 62, 62, 62, 62, 62, 62, 63, 63, 63, 63, 63, 63, 63, 63, 64, 64, 64, 64, 64…
$ LATITUD_MINUTOS <int> 52, 52, 52, 52, 52, 52, 52, 52, 14, 14, 14, 14, 14, 14, 14, 14, 24, 24, 24, 24, 24, 24, 24, 24, 14, 14, 14, 14, 14…
$ LONGITUD_GRADOS <int> 34, 34, 34, 34, 34, 34, 34, 34, 58, 58, 58, 58, 58, 58, 58, 58, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56…
$ LONGITUD_MINUTOS <int> 34, 34, 34, 34, 34, 34, 34, 34, 38, 38, 38, 38, 38, 38, 38, 38, 59, 59, 59, 59, 59, 59, 59, 59, 43, 43, 43, 43, 43…
$ ALTURA <int> 256, 256, 256, 256, 256, 256, 256, 256, 11, 11, 11, 11, 11, 11, 11, 11, 24, 24, 24, 24, 24, 24, 24, 24, 198, 198, …
$ NRO <int> 89034, 89034, 89034, 89034, 89034, 89034, 89034, 89034, 89053, 89053, 89053, 89053, 89053, 89053, 89053, 89053, 88…
$ NroOACI <chr> "SAYB", "SAYB", "SAYB", "SAYB", "SAYB", "SAYB", "SAYB", "SAYB", "SAYJ", "SAYJ", "SAYJ", "SAYJ", "SAYJ", "SAYJ", "S…
$ LATITUD <dbl> -77.86667, -77.86667, -77.86667, -77.86667, -77.86667, -77.86667, -77.86667, -77.86667, -62.23333, -62.23333, -62.…
$ LONGITUD <dbl> -34.56667, -34.56667, -34.56667, -34.56667, -34.56667, -34.56667, -34.56667, -34.56667, -58.63333, -58.63333, -58.…
$ FECHA <int> 20042021, 20042021, 20042021, 20042021, 20042021, 20042021, 20042021, 20042021, 20042021, 20042021, 20042021, 2004…
$ HORA <int> 0, 3, 6, 9, 12, 15, 18, 21, 0, 3, 6, 9, 12, 15, 18, 21, 0, 3, 6, 9, 12, 15, 18, 21, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, …
$ TEMP <dbl> -20.0, -16.2, -16.9, -13.0, -15.1, -16.6, -15.3, -16.8, -0.1, -0.5, -0.9, -0.5, 0.3, 1.1, 1.3, 2.5, 3.9, 3.1, 3.5,…
$ HUM <int> 70, 72, 80, 73, 86, 66, 77, 85, 88, 82, 77, 75, 84, 92, 99, 91, 43, 49, 67, 81, 67, 62, 77, 83, 74, 78, 75, 75, 74…
$ PNM <dbl> 975.7, 970.4, 966.7, 964.2, 961.4, 963.0, 964.8, 967.1, 992.9, 993.0, 992.9, 994.6, 995.3, 995.6, 994.9, 992.3, 98…
$ DD <int> 140, 140, 110, 140, 160, 160, 90, 110, 230, 230, 270, 290, 290, 320, 320, 320, 160, 200, 50, 200, 110, 250, 250, 2…
$ FF <dbl> 13, 30, 50, 52, 74, 61, 44, 13, 15, 15, 26, 28, 41, 37, 50, 74, 20, 7, 13, 7, 11, 46, 56, 74, 19, 20, 19, 20, 15, …
Eliminamos las observaciones nulas en el dataset
data = na.omit(data)
corr <- cor(data[, c(7,13,14,15,16,17,18)], use = "complete.obs")
ggcorrplot(corr, type = "lower", outline.col = "black",
lab=TRUE,
ggtheme = ggplot2::theme_gray,
colors = c("#6D9EC1", "white", "#E46726"))
Agregamos los datos por nombre calculando el promedio y desvio del viento
data_agg = data %>%
group_by(NOMBRE) %>%
summarise(MEAN_VIENTO_KMH = mean(FF), SD_VIENTO_KMH = sd(FF), LONGITUD = unique(LONGITUD), LATITUD = unique(LATITUD), .groups = "keep")
Dibujamos en el mapa de argentina los puntos que tenemos en el dataset agrupado. Para eso cargamos el dataset de departamentos que tiene info de toda la argentina
library(ggplot2)
departamentos <- st_read("data_departamentos/Codgeo_Pais_x_dpto_con_datos/pxdptodatosok.shp")
Reading layer `pxdptodatosok' from data source `/Users/antonellaschiavoni/Antonella/Maestria/EstadisticaEspacial/TPFinal/Datos/trabajo_final/smn/estadistica_espacial/data_departamentos/Codgeo_Pais_x_dpto_con_datos/pxdptodatosok.shp' using driver `ESRI Shapefile'
Simple feature collection with 527 features and 10 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -74.02985 ymin: -90 xmax: -25.02314 ymax: -21.74506
Geodetic CRS: WGS 84
#verificamos la proyeccion de departamentos
st_crs(departamentos)
Coordinate Reference System:
User input: WGS 84
wkt:
GEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["latitude",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["longitude",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]]
#como tiene otra proyeccion la pasamos a mercator
departamentos <- st_transform(departamentos, crs = 4326)
#ggplot() +
# geom_sf(data = departamentos) +
# geom_point(data = data_agg, aes(x = x, y = y), colour = "red", size = 1) +
# theme_minimal() +
# ggsave('plot.png',height = 25 , width = 30)
Convertimos data_agg a data.frame
df_data_agg = data.frame(data_agg)
Transformamos df_data_agg en un archivo geográfico utilizando el código de proyección mercator
data_sf = st_as_sf(df_data_agg, coords = c("LONGITUD", "LATITUD"), crs = 4326)
Validamos la clase del nuevo dataframe
class(data_sf)
[1] "sf" "data.frame"
##TODO: graficar el mapa de densidades con los puntos de las estaciones en data.
library(leaflet)
leaflet() %>%
addTiles() %>%
addCircles(data = data_sf, weight = 3)
NA
Analisis univariado
describe(data_agg$MEAN_VIENTO_KMH)
hist(data_agg$MEAN_VIENTO_KMH)
boxplot(data_agg$MEAN_VIENTO_KMH)
hist(data_agg$SD_VIENTO_KMH)
boxplot(data_agg$SD_VIENTO_KMH)
hist(data_agg$MEAN_VIENTO_KMH)
boxplot(data_agg$MEAN_VIENTO_KMH)
qqnorm(data_agg$MEAN_VIENTO_KMH)
qqline(data_agg$MEAN_VIENTO_KMH, col=2)
hacemos el test de normalidad
shapiro.test(data_agg$MEAN_VIENTO_KMH)
Shapiro-Wilk normality test
data: data_agg$MEAN_VIENTO_KMH
W = 0.96947, p-value = 0.02215
El p-value obtenido es menor a 0.05, lo cual indica que los datos que tenemos no son normales
coordinates(data_agg) = ~LATITUD+LONGITUD
class(data_agg)
[1] "SpatialPointsDataFrame"
attr(,"package")
[1] "sp"
plot(data_agg)
vt_exp = fit.variogram(v, vgm(125, "Exp", 30, 5))
vt_exp
plot(v , vt_exp)
vt_mat = fit.variogram(v, vgm(125, "Mat", 30, 5))
plot(v , vt_sph)
vt_exc = fit.variogram(v, vgm(125, "Exc", 30, 5))
plot(v , vt_sph)
attr(vt_exp, 'SSErr')
[1] 0.0126554
attr(vt_mat, 'SSErr')
[1] 0.0126554
attr(vt_exc, 'SSErr')
[1] 0.01653322